
data_1 %<>%
  mutate(income_recode2 = case_when(q26_survey_income_recode %in% c(1, 2) ~ 1,
                                    q26_survey_income_recode == 3 ~ 2,
                                    q26_survey_income_recode == 4 ~ 3,
                                    q26_survey_income_recode == 5 ~ 4,
                                    q26_survey_income_recode == 6 ~ 5,
                                    q26_survey_income_recode %in% c(7, 8) ~ 6,
                                    TRUE ~ NA_real_),
         q22_survey_working_recode = case_when(q22_survey_working == "Ya" ~ 1,
                                               q22_survey_working == "Tidak" ~ 0,
                                               TRUE ~ NA_real_),
         q23_survey_look_work_recode = case_when(q23_survey_look_work == "Ya" ~ 1,
                                                 q23_survey_look_work == "Tidak" ~ 0,
                                                 TRUE ~ NA_real_)
  )





test_mod1a <- lm(income_recode2 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod2a <- lm(q22_survey_working_recode ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod3a <- lm(q27_survey_job_satis ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod4a <- lm(income_recode2 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5), got_job == 0)
test_mod5a <- lm(q22_survey_working_recode ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5), got_job == 0)
test_mod6a <- lm(q27_survey_job_satis ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5), got_job == 0)


observations <- c(nobs(test_mod1a), nobs(test_mod4a), 
                  nobs(test_mod2a), nobs(test_mod5a), 
                  nobs(test_mod3a), nobs(test_mod6a))

test_mod1a <- coeftest(test_mod1a, vcov=vcovHC(test_mod1a,type="HC0"))
test_mod2a <- coeftest(test_mod2a, vcov=vcovHC(test_mod2a,type="HC0"))
test_mod3a <- coeftest(test_mod3a, vcov=vcovHC(test_mod3a,type="HC0"))
test_mod4a <- coeftest(test_mod4a, vcov=vcovHC(test_mod4a,type="HC0"))
test_mod5a <- coeftest(test_mod5a, vcov=vcovHC(test_mod5a,type="HC0"))
test_mod6a <- coeftest(test_mod6a, vcov=vcovHC(test_mod6a,type="HC0"))



#Javanese advantage table

table <- list(test_mod1a, test_mod4a ,test_mod2a,test_mod5a, test_mod3a, test_mod6a)

note_text <- paste(" Beta coefficients from OLS regression with interaction terms. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Basic Competence Examination (SKD) Failure On Work Outcomes",
                  label = 'tab:testing_effect_work',
                  model.names = F,
                  model.numbers = T,
                  digits = 2,
                  column.separate = c(2, 2, 2),
                  column.labels = c("Income (m, IDR)", "Employed", "Job Satisfaction (1-4)"),
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(c("Subset", "----", "Non-CS", "---", "Non-CS", "---", "Non-CS"),
                                   c("Observations", observations),
                                   c('Bandwidth', rep(c('1\\%'), 6))),
                  covariate.labels = c("Failed Test"),
                  star.cutoffs = c(0.05, 0.01, 0.001),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 28)], note_text, './_4_outputs/tables/table_a15.tex', .8)


